Density, structure and dynamics of water: the effect of van der Waals interactions 
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It is known that ab initio molecular dynamics (AIMD) simulations of liquid water, based on the 
generalized gradient approximation (GGA) to density functional theory (DFT), yield structural and 
diflusive properties in reasonable agreement with experiment only if artificially high temperatures 
are used in the simulations. The equilibrium density, at normal conditions, of DFT water has been 
recently shown by Schmidt et al. [J. Phys. chem. B, 113, 11959 (2009)] to be underestimated by 
different GGA functionals for exchange and correlation, and corrected by the addition of interatomic 
pair potentials to describe van der Waals (vdW) interactions. In this contribution we present a DFT- 
AIMD study of liquid water using several GGA functionals as well as the van der Waals density 
functional (vdW-DF) of Dion et al. [Phys. Rev. Lett. 92,246401(2004)]. As expected, we find 
that the density of water is grossly underestimated by GGA functionals. When a vdW-DF is used, 
the density improves drastically and the experimental diffusivity is reproduced without the need of 
thermal corrections. We analyze the origin of the density differences between all the functionals. 
We show that the vdW-DF increases the population of non-H-bonded interstitial sites, at distances 
between the first and second coordination shells. However, it excessively weakens the H-bond 
network, collapsing the second coordination shell. This structural problem is partially associated 
to the choice of GGA exchange in the vdW-DF. We show that a different choice for the exchange 
functional is enough to achieve an overall improvement both in structure and diffusivity. 
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I. INTRODUCTION 

Studies with density functional theory (DFT)-based 
ab initio molecular dynamics (AIMD) simulations of liq- 
uid water-ir— seriously differ from experimental results 
both in diffusive and structural properties. This is true 
even though most studies have been performed at the 
experimental density, what implies theoretical pressures 
as high as 1 GPa4 The origin of the discrepancy with 
experiments2r2i is still unclear. One of the main con- 
clusions of these studies^ was that the AIMD results for 
radial distribution functions (RDF) and self-diffusivity at 
a given temperature compare well with the experimental 
results at a temperature 20% lower. Therefore, a prac- 
tical solution to the problem is to perform simulations 
at temperatures 20% higher than the desired reference 
temperaturCjiS as proposed in ref y. However, although 
this ad hoc solutions may be convenient to study other 
properties of liquid water^i^ of molecules in solution^i^ 
and of wet interfacesr^ they do not help explain the ori- 
gin of the differences between experiments and AIMD 
simulations. 

There are obvious limitations in the AIMD description 
of liquid water that could account for these differences, 
including the inability of present generalized-gradient ap- 
proximation (GGA) density functionals to describe dis- 
persion interactions, or the complete neglect of quantum 
fluctuations in the classical treatment of nuclear dynam- 
ics. The latter question has already been addressed: 
Morrone and Car^^' have found that when nuclear quan- 
tum effects are accounted for, the simulated liquid be- 
comes less structured and more diffusive, and the agree- 



ment between AIMD simulations and experiments con- 
siderably improves. On the other hand, this result is 
at odds with a previous similar study in which Chen et 
ati^ found that the quantum vibrations of the hydrogen 
atoms would indeed increase the strength of the hydro- 
gen bonds in liquid water. Recently, Habershon et a/— 
identified the competing contributions from intermolec- 
ular and intramolecular quantum fluctuations to cancel 
out and thus to have a small net effect on the diffusion 
coefficient. 

Computational studies of liquid water differ widely in 
their methodology: e.g. they might employ either Born- 
Oppenheimer or Car-Parrinelloi^ dynamics, they can use 
different types of basis setSfii^i^ii^, and nuclear quantum 
effects may or may not be accounted for. Nonetheless, 
most are performed at a fixed volume, set approximately 
to the experimental density of water at room temper- 
ature, ~1.0 g/cm'^. This approach would be correct if 
the calculated pressure also reproduced the experimen- 
tal pressure at normal conditions (1 atm). However, as 
pointed out in Refs. |4||20|. even though pressure fluctua- 
tions are very large, due to the small size of the simulated 
system, the observed average pressure is much larger than 
the experimental one. Therefore, both structural and dif- 
fusive features correspond to a region of the temperature- 
pressure phase diagram which differs form the experi- 
mental values they are compared to32. In a word, the 
simulations are performed at too high temperatures and 
pressures. 

Schwegler et al^ have already presented a study of 
DFT water under pressure, but the equilibrium density of 
GGA water was not explored. In Ref.[20i, the discrepancy 



between theory and experiments of the pressure-density 
phase diagram was aheady addressed using Monte Carlo 
simulations. A recent study by Schmidt et aL- of liq- 
uid water under isothermal and isobaric condition, us- 
ing both GGA and so-called DFT+dispersion (DFT- 
D) method has analyzed this question in detail. They 
have shown that two of the most commonly used GGAs, 
BLYP^^'^'^ and PBE^* have an equilibrium density of 
0.75 and 0.88 g/cm^ respectively, 25% and 12% below 
the experimental value. When vdW interactions are 
added, using the interatomic pair potentials proposed by 
Grimme)^ the experimental density is recovered. 

In this study, we partly confirm and extend their re- 
sults by using AIMD but, more importantly, we provide 
insights into the origin of the vdW effects: by compar- 
ing results with the GGA and with the first-principles 
vdW density functional (vdW-DF) of Dion et al,^^ and 
by separating H-bonded and non-H-bonded interactions, 
we are able to explain the differences in density and local 
structure between GGA water and vdW-DF water. 



II. METHODOLOGY 

The simulations are performed using the self-consistent 
Kohn-Sham approach^^ to DFT— within the generalized- 
gradient approximation to exchange and correlation 
(XC). We choose two commonly used GGA functionals: 
PBE^^ and rcvPBE.— The choice of these functionals is 
motivated by our interest on studying the effect of dis- 
persion interactions using the vdW-DF<^ The original 
vdW-DF, that we coin DRSLL after its authors, ^^ uses 
revPBE exchange combined with the non-local vdW cor- 
relation. Since this choice was somewhat arbitrary, we 
have additionally substituted it by PBE exchange, label- 
ing this functional DRSLL-PBE. Besides addressing the 
effect of the exchange part in vdW-DF, this also pro- 
vides a second reference point on the effect of non local 
correlations on GGA functionals, i. e. DRSLL vs revPBE 
and DRSLL-PBE vs PBE. We have used the implemen- 
tation of vdW-DF by Roman-Perez and Soler'^" within 
the SIESTA program. 



result is presented in Table HI We provide the details of 
these two sets in the Supplementary Information. 

Our results show that the main difference between the 
DZP basis set we have used in this study and the well 
converged TZP basis is that the smaller basis produces 
a somewhat less structured liquid. The diffusivity is ap- 
proximately the same for the two basis sets. The first 
peak in the 0-0 radial distribution function (RDF) is 
almost identical, the first minimum is 13% higher. The 
comparison of the radial distribution function between 
the TZP basis and the DZP basis is shown in Fig.[T] The 
error of the DZP basis on the location of the first maxi- 
mum in the 0-0 RDF is less than 0.6%, indicating that 
the basis is very good in the description of the H-bond 
geometry, as already shown in Ref. 4. 

In order to evaluate total energy errors associated to 
our choice of basis set, we have selected 10 snapshots, sep- 
arated by 2 ps time intervals, from our simulation density 
1.00 g/cm^. The difference between the TZP and DZP 
energies is an almost constant value of 0.12 eV/molecule, 
with a standard deviation of only 7 meV, which can be 
neglected even when we are studying the small van der 
Waals effects. Also, in Table HI the details of the two 
simulations with DZP and TZP are compared. The en- 
ergy drift is on the order of 10~^ a.u./atom/ps, which 
accounts for a change of only a few Kelvin during 20 
ps, well below the statistical fluctuation of temperature 
from the MD simulation. We have estimated that basis- 
set-induced errors in the equilibrium pressure are on the 
same order of the error bar associated to the choice of 
different initial configurations, which is approximately 1 
kbar. Both errors are much smaller than the statistical 
fluctuations of pressure during the MD simulation, which 
is on the order of 3-4 kbar. Overall, we can conclude 
that AIMD of water are not totally insensitive to the ba- 
sis set choice,— li^ shorter basis sets providing slightly less 
structured liquids in general. However, the results and 
conclusions presented in this study are well beyond the 
uncertainties due to these limitations. 



A. Basis set 

Core electrons were replaced by norm- 
conserving pseudopotentials^ in their fully nonlocal 
representation.'^^ Numerical atomic orbitals of finite 
support were used as basis set, and the calculation of 
the self-consistent Hamiltonian and overlap matrices 
was done using the siesta method.— i2i In this study, 
a double-C polarized (DZP) basis set was used, which 
has been variationally optimized following the method 
proposed in Refs. [35ll36l The validation of the method, 
pseudopotentials and basis set can be found in Ref. 0. 
In addition, a triple-C polarized (TZP) basis set has also 
been tested in this study. The summary of the basis test 



TABLE I: Summary of AIMD simulations results, for liquid 
water at 1.0 g/cm"^ and 340 K, using the PBE functional. 
We compare results with a DZP basis using a confining pres- 
sure of 0.2 GPa with a run of 20 ps, and a TZP basis using 
a confining pressure of 0.2 GPa with a run of 15 ps. Tem- 
perature T (K), pressure P (kbar), energy drift during the 
simulation Edrift (10~^ a.u./atom/ps), self diffusion coeffi- 
cient D (10"^cm^/s), height of first peak Qoo' , height of first 
minimum go(3 and position of first maximum r^o^ (A) of 
the O-O radial distribution function are compared. 
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FIG. 1: Comparison of the O-O radial distribution function 
(RDF) obtained with the DZP(sohd black) and TZP(dashed, 
red) basis sets (see text) at 1.0 g/cm'^. 



B. Molecular dynamics simulations 

All the results are based on AIMD simulations of 
64 molecules of heavy water. Classical molecular dy- 
namics simulations, where the electronic potential seen 
by the nuclei is replaced by an empirical force field, 
were performed before AIMD equilibration in order to 
prepare reasonably equilibrated initial configurations. 
All empirical-potential-based simulations were performed 
with the TIP4P force field^ as implemented within 
the GROMACS MD package^*-?^ under constant vol- 
ume and temperature conditions, with a Nose-Hoover 
thermostat^Siii at 300 K, along a 1 ns equilibration tra- 
jectory. AIMD simulations were started from these pre- 
equilibrated systems imposing an additional 3 ps equi- 
libration by means of temperature annealing (velocity 
rescaling).— The actual production runs of the AIMD 
simulations (20 ps each) were accomplished by constant- 
energy Verlet's integration, given our interest in dynam- 
ical properties. The time step in all simulations, (includ- 
ing empirical force field MD and AIMD) is 0.5 fs. 

The simulations are performed under constant volume 
i. e. fixed cell size and shape, under periodic boundary 
conditions. Rather than performing a constant-pressure 
simulation to establish the theoretical equilibrium den- 
sity under ambient conditions, we chose to perform a se- 
ries of simulations by changing the volume of the unit cell 
while keeping the number of water molecules fixed. The 
reduced system size produces large pressure and tempera- 
ture fluctuations and requires larger statistics for reliable 
comparisons between different simulations. Therefore, 
we compute three different trajectories for each density, 
starting from different equilibrations. In total, over 500 
ps of AIMD simulations were produced in this study. 

The temperature of the simulations using GGA func- 
tionals were set higher than 300 K due to the reason 
mentioned in the introduction. With the PBE functional 



we use ~360 K^ and with revPBE we use ^^330 K4^ 



C. System size effects 

System size effects on structural and diffusive prop- 
erties of simulated liquid water have been addressed 
before. Most studies^'^^ used empirical potentials and 
found that size effects seem not to be a problem for those 
properties. In AIMD simulations, size effects have been 
shown to be rather small for structural properties )^i^ 
but much larger for dynamical properties^^i^meaning 
that finite size scaling is needed to obtain the infinite-size 
self-diffusivity coefficient^^ With such finite size scaling, 
Kuhne et al found^ an improved self-diffusivity but still 
smaller than experiment. 

We have explored the influence of size-dependent 
boundary conditions on both structural and dynamical 
properties of the system. The dependence of average 
pressure on the system size has also been studied in de- 
tail. These analyses were performed by means of classical 
molecular dynamics simulations. It seems reasonable to 
assume that size effects are robust enough to the change 
of interaction potential and that results obtained from 
force-field-based simulations are adequate to estimate the 
error due to the reduced size of the system in our AIMD 
studies. 

We have performed molecular dynamics simulations 
of five different system sizes (32, 64, 128, 512 and 1024 
molecules) at four different densities (0.95, 1.0, 1.1, and 
1.2 g/cm^) with the empirical force field TIP4P.'^'^ These 
simulations were equilibrated at 300 K along 100 ps 
runs using a Nose-Hoover thermostat. Subsequently, 
they were allowed to continue for 1 ns, using a Verlet 
integrator)^ 
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FIG. 2; Convergence of average pressure with system size. 
AP is the pressure difference between the system of A*' 
molecules and that of 1024 molecules, at 300 K and four differ- 
ent densities. All the simulations in this figure were performed 
using the TIP4P empirical force field. 



The convergence of average pressure of TIP4P simu- 
lations with system size is shown in Fig. [2] Simulations 
with 32 water molecules result in a large discrepancy with 
the converged value. With 64 molecules the pressure is 
not yet fully converged but the error is much smaller. 
Considering the need for long simulation times and the 
large number of simulations required for this study, we 
chose to limit our systems to 64 water molecules to study 
the density-pressure dependence. 



III. SIMULATION RESULTS 
A. Equilibrium density for GGA functionals 

Table |TT] summarizes the parameters used, and the av- 
erage values obtained, in all our AIMD simulations, span- 
ning five functionals and a wide range of pressures and 
densities. As explained before, the simulation tempera- 
tures were chosen, for each functional, to achieve a rea- 
sonable comparison with the experimental difFusivities 
and radial distribution functionsi^i^ 

Pressure-density curves obtained with three GGA 
functionals (BLYP, PBE and revPBE) and two vdW- 
DFs (DRSLL and DRSLL-PBE) are compared in Fig. |3l 
Because of the reasons given in the introduction, each 
functional is simulated at a different temperature (see 
Table HIl) . but wc would like to stress that this has only 
a minor effect on the resulting average pressures 4ii^ To 
illustrate this, we show in Fig. [3] the experimental densi- 
ties at 300 K and 360 K for P = 1 atm. Their difference 
is much smaller than that between different functionals. 
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densities considerably lower than the experimental value 
(p — 0.997 g/cm'^ at P = 1 atm). The equilibrium 
density of PBE is 0.85-0.90 g/cm^, close to the value 
presented in Refs. 8,21 and around 12% lower than ex- 
periments. The BLYP equilibrium density of 0.76-0.85 
g/cm'^, also within the range provided by Schmidt et a& 
and 19% lower than experiments. The revPBE density 
is 0.63-0.75 g/cm'^, or 31% less dense than experiments. 

The large differences in the calculated densities of sim- 
ilar GGAs are surprising and require some analysis. Like 
many liquids, water maintains much of the short-range 
order of its solid. In ice Ih, a rigid framework of hydro- 
gen bonds forces a tetrahedral coordination and a rel- 
atively open structure with large voids. Four of these 
interstitial voids surround each molecule, along orienta- 
tions opposite to those of its H bonds (we will call them 
'anti-tetrahedral' orientations). In the liquid, entropy im- 
plies that part of the H bonds are broken and part of the 
interstitial voids are occupied. Thus, we can imagine 
three different properties of a functional or force field, 
that will determine the density of liquid water: i) the 
length of the hydrogen bonds, causing a possible dila- 
tion of the H-bond framework; ii) the strength of these 
bonds, that determines the average H-bond coordination; 
and in) the potential energy of molecules at interstitial 
positions, that determines their average occupation. 

With a rigid H-bond framework, the density would be 
proportional to doc where doo is the average 0-0 dis- 
tance. The change in doo between PBE and revPBE (see 
Table [n| is ~3%, which would translate into a difference 
in density of just ~10%, far from the ^25 % observed. 

However the density also depends on the number of 
H bonds, which increases with their strength. Fig. ^h) 
shows that, despite the similarity of the two functionals, 
the H bond energy is ~50% larger with PBE than with 
revPBE. As a consequence, the first coordination peak 
in the 0-0 radial distribution function (RDF) goo{'>'), 
plotted in Fig. 21 is ^^30% higher with PBE, showing 
that the H-bond network is better preserved with this 
functional. All put together, the difference in densities 
within the first coordination shell, calculated as ,"^°°°iJ"L 

(with values from Table Hi]) , where Ncoor is the average 
coordination number and t-qo" is the radius of the first 
minimum in gooif): is larger than 30%. 
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B. GGA functionals: effect of pressure on the RDF 



FIG. 3: Pressure-density curves obtained in AIMD simula- 
tions with different GGA and vdW XC functionals. The er- 
ror bars are standard deviations of average pressures in three 
different runs with the same density, XC functional, and tem- 
perature. The simulation temperatures can be found in Table 
im The colored boxes show the estimated range of equilibrium 
densities at P = 1 atm ~ 1 bar. The arrows indicate the ex- 
perimental densities at 300 K and 360 K and 1 atm.— 



All GGA functionals result in theoretical equilibrium 



Fig. |4] compares the 0-0 RDF of liquid water at 0.95 
g/cm'^ for the PBE and revPBE functionals. The temper- 
ature difference between the the two simulations is under 
15 K, and should have a negligible effect on the RDFs^i 
compared with that resulting from pressure differences. 
We compare our results with two different experimental 
data sets^i^, both obtained from Ref. IM The revPBE 
RDF matches the experiments rather well. Although the 
onset and position of the first peak are clearly shifted to 
the right, there is a large uncertainty in the experimental 



TABLE II: Simulations parameters used, and average values obtained in this work: mass density p, exchange-correlation 
functional, average temperature Tavg, diffusion coefficient D, position tqq' of first maximum in goo{r), position roo" of first 
minimum in goo{r), average coordination number, and average number of hydrogen bonds A^'H-bond- For each calculation, 
quantities are averaged over 20 ps. At ambient conditions, the experimental self-diffusion coefficient is 2.2 x 10~^ cm^/s for 
H2O and 1.8xl0~ cm /s for D2O-. The experimental coordination number is 4.7^^ 



P (g/cm'') 


Functional 


Ta., (K) 


D (lO^^'cm^s) 


ro5^ (A) 


'■^00 (A) 


Coord. 


A^H-bonds 


0.65 


BLYP 


361 


2.82 
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0.85 
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352 


1.33 


2.77 


3.41 


4.28 


3.67 
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354 


1.16 
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3.36 


4.31 


3.69 


1.00 


PBE 
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1.18 


2.75 


3.28 


4.23 


3.78 


0.65 
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4.35 
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3.62 


3.58 
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0.75 
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0.95 
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DRSLL 
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3.63 


1.00 
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3.46 


4.90 
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2.83 
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"As there is no sharply identifiable first minimum in goo{i') at 
1.05 g/cm'^, we can not provide the accurate value here. 



data for this peak (see Ref. l50|). However, at 0.95 g/cm'^, 
revPBE water is at a very large pressure. When this 
pressure is released, the density drops to '^0.65 g/cm'^ 
and the agreement with the experimental RDF is totally 
lost (see Supplementary Information). Both the first and 
second peaks move further to the right. This effect is also 
observed in PBE. It should be noted that, for PBE, the 
height of the first minimum remains well below the ex- 
perimental value and is barely modified, either by tem- 
perature or pressure effects (as seen in supplementary 
information) . 



C. Van der Waals Interactions in Water 

Van der Waals (vdW) dispersion interactions, due to 
non local electron correlations, are not treated properly 
by the local density^ (LDA) and GGAs^ii^iS function- 
als, in which electron correlations are treated as local or 
semi-local effects. It is known that in water, due to the 
high polarizability of oxygen^^"— vdW interactions have 
a significant contribution to the binding. The vdW at- 
traction contributes to strengthening both H-bonds and 
non-H-bond interactions' and it increases the overall co- 
hesive energy in the liquid. Schmidt et al have studied 
the vdW effect on the density of water with the PBE-D 
method,— which includes an interatomic pair potential 
correction added to the PBE functionali^ They showed 
that the density of PBE-D water is very close to the 
experimental value, and the resulting liquid is also struc- 
turally closer to experiments. The dynamical properties 
of their liquid were not accessible due to the use of a tem- 
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FIG. 4: 0-0 radial distribution functions (RDF) of liquid 
water for AIMDs in this work, at 0.95 g/cm^ with PBE at 
360 K (solid black line) with revPBE at 340 K (solid red 
line). These temperatures are chosen^ so that the structure 
and diffusivity compare best with those of experiments at 300 
K. Expl (dotted line) and Exp2 (dashed line) data are the 
experimental data in Refs. [43 and[50| respectively. 



perature and pressure thermostat in their simulations. 

Furthermore, classical force fields, that represent elec- 
tronic dispersion correlations with — l/r^ interatomic 
pair potentials, can be inaccurate because these corre- 
lations are not atom-centered in principle and they de- 
pend on the instantaneous environment of each atom and 
molecule. 



Since the current DFT description of liquid water over- 
estimates intermolecular binding, it is not clear whether 
vdW interactions will reduce the difference with exper- 
iments of various magnitudes. However, the attractive 
vdW effects could be especially relevant to increase the 
density of the simulated liquid, bringing it closer to ex- 
periment, as Schmidt et al have seen using the PBE-D 
method. 

As seen in Fig. |3l the equilibrium density computed 
with the DRSLL vdW-DF is 1.02 g/cm^, only 2% above 
the experimental value. This is much better than any of 
the commonly used GGA functionals. 

The 0-0 RDF of DRSLL water show large disagree- 
ments with experiments, as shown in Fig. [5l 0-H and 
H-H RDFs are presented in the supplementary infor- 
mation. The average first neighbor distance, is longer 
than the experiments. The second coordination peak in 
<?oo('') 'collapses' significantly with this functional, and 
an anomalous hump at r = 3.8 A appears. 

Differences in the liquid H-bond interactions can 
be characterized by relaxing typical liquid-phase 
structuresi^ Table Hill shows several characteristic mag- 
nitudes calculated after relaxing five pentamer clusters, 
randomly chosen from our AIMD,^ with different func- 
tionals (but with the same initial geometries for all func- 
tionals). The table shows that DRSLL increases the 
H-bond binding energy of revPBE but it also yields a 
longer H-bond, producing the shift of the first coordina- 
tion peak to longer distances. This overestimation of 
DRSLL binding distances has been observed in many 
other systems <^i^ 



TABLE III: Comparison of average H-bond energy Ehb, 
and distance (Ihb, first neighbor O-O distance doc, and in- 
tramolecular 0-H bond length doH obtained by relaxing the 
same five pentamer clusters with different XC functionals. 
Each cluster contains a randomly-chosen central molecule plus 
its first neighbors. 

Functional Ehb (eV) duB (A) dpo (A) dpH (Xf 



BLYP 


0.222 


1.921 


2.881 


0.978 


FEE 


0.262 


1.872 


2.809 


0.979 


revPBE 


0.199 


2.022 


2.974 


0.974 


DRSLL 


0.225 


2.022 


2.980 


0.974 


DRSLL-PBE 


0.290 


1.918 


2.833 


0.976 



D. Effect of PBE exchange in vdW-DF 

As mentioned before, the exchange energy in 
the DRSLL functional is taken from the revPBE 
formulation.'^^ This choice reduces, relative to PBE, the 
spurious exchange-induced binding in simple systems like 
noble gas dimersi^ However, revPBE underestimates the 
H-bond interaction, what translates to larger 0-0 dis- 
tances and a clear under-structuring of the RDF, com- 
pared to experiments. Therefore, it is worth studying the 
effect of combining vdW correlations, as described by the 
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FIG. 5: Comparison of the O-O RDFs obtained with the 
DRSLL (black solid line) and DRSLL-PBE (red solid fine) 
functionals, and experimental results of Refs. 49 (dotted line) 
andlSOl (dashed line). Both simulations and experiments are 
at 1.0 g/cm^ and ~300 K. 



DRSLL functional, with a stronger H-bond interaction, 
as obtained using PBE exchange. Thus, we have replaced 
the revPBE exchange in DRSLL by PBE exchange and 
we will refer to this functional as DRSLL-PBE. Recent 
studies^°-^^ have shown that minor modifications of the 
exchange enhancement factor, in the same line as that 
used here, can improve geometries and binding energies 
of small molecular systems. However in this study we 
chose the original DRSLL and DRSLL-PBE to allow for 
comparisons with their GGA counterparts. 

We have performed simulations with DRSLL-PBE for 
1.00, 1.05, 1.10 and 1.20 g/cm^ at 300 K. The data points 
in Fig. [3] for this functional represent a single simula- 
tion each, and we used the same error bars obtained for 
DRSLL. At 1.00 g/cm^ average pressures are ^10 kbar, 
indicating that the strong binding of PBE exchange, 
added to that induced by vdW correlations, has a net 
effect of over correcting the positive pressure of GGAs. 
We find an equilibrium density at 1 atm of 1.13 g/cm'^ In- 
deed, the pressure reduction achieved by DRSLL versus 
revPBE is 13 kbar, close to that of DRSLL-PBE versus 
PBE, 14 kbar. 

As shown in Fig. [3 DRSLL-PBE recovers an 0-0 RDF 
similar to experiments, without an artificial increase in 
the simulation temperature. The height of the first min- 
imum in 5oo(^) is very close to the experimental one, 
and the anomalous hump displayed by DRSLL is almost 
completely eliminated. The fact that the long range tail 
of goo {r) matches closely the experimental one may be 
indicative of a correct characterization of the long range 
structure of water, obtained when vdW interactions are 
included. We provide a comparison of the structure fac- 
tors calculated for PBE and DRSLL-PBE in the supple- 
mentary information. 



E. Self Diffusion Properties 

In Fig. [6l the temperature dependence of diffusivity, 
computed from all our AIMD simulations, is compared 
with experimental values at similar and lower tempera- 
tures. We also show values corrected for finite size effects, 
as proposed in Ref. [46!. Confirming previous results,- we 
find that the temperature of AIMD with PBE and BLYP 
must be 16-20% higher than the experimental one, to 
give a similar diffusivity. With revPBE, the temperature 
needs to be only 6-9% higher. 
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FIG. 6: Self-diffusion coefficient vs temperature obtained in 
AIMD simulations with various XC functionals compared to 
H2O and D2O experimental values.— Values corrected for fi- 
nite size effects^ are shown as soUd symbols, and observed 
values as open symbols. Note that our simulations are per- 
formed for D2O. 



On the other hand, both DRSLL and DRSLL-PBE 
nearly reproduce the experimental diffusivity at room 
temperature, without any temperature rescaling. This 
is not unexpected, given the strong link between struc- 
ture and dynamical properties of liquid water j^ Thus, 
one of our most important conclusions is that DRSLL- 
PBE is a very good XC functional to describe liquid wa- 
ter: PBE requires scaling up the temperature to reach an 
agreement with experiments. DRSLL replicates diffusiv- 
ity well, but it substantially fails to describe the RDFs. 
DRSLL-PBE corrects the structure while maintaining a 
very good diffusivity. 



IV. ANALYSIS AND DISCUSSION 
A. Intermolecular interactions 

Since the vdW attraction of vdW-DFs does not re- 
duce the H-bond distances (actually they increase, rela- 
tive to their GGA counterparts), what causes the large 
increase in equilibrium density? One possibility is that 
the stronger H-bonds, shown in Fig. ^h), increase the 
average first-neighbor coordination. However, the height 



of the first peak in gooif), as shown in Figs. |4] and [5l ac- 
tually decreases from the GGAs to the vdW-DFs. There- 
fore, we are left only with the third possibility mentioned 
before, i. e. an increase in the occupation of the 'intersti- 
tial', non-H-bonded sites. 

In order to characterize the vdW binding between wa- 
ter molecules, we have calculated the interaction en- 
ergy between two water molecules oriented as shown in 
Fig. [7l^a), thus avoiding the hydrogen bond interaction. 
By comparison, we also use PBE and revPBE. Fig. El^a) 
clearly shows that DRSLL exhibits a minimum in the po- 
tential at ^3.7 A with a binding energy of 10 meV, not 
shown by the GGAs. Although this is many times times 
weaker than the H-bond, it will have an important effect 
in increasing the occupation of the interstitial sites, which 
are roughly at that distance. It is worth noting that the 
position of this potential energy minimum is close to the 
first minimum of 500(^)7 obtained with the GGAs. This 
means that this new 'vdW bond', and its effect on in- 
creasing the occupation of interstitial sites, may account 
for the unusual hump appearing in the goo (i^) of DRSLL 
at this distance, as seen in Fig. [5j The potential energy 
minimum in DRSLL-PBE is even deeper, 25 meV, and 
shifted to shorter distances (3.4 A). This is also the origin 
of the increase of the height of the first gooir) minimum 
from PBE to DRSLL-PBE. 

We have also computed the H-bond potential energy 
curve between two water molecules in a H-bound config- 
uration. The results are shown in Fig. [7l^b). The depth 
of the potential is lower that the optimal H-bond inter- 
action energy because the geometry of the molecules was 
not optimized. Still the results show the same tendency 
observed in the RDF and in Table Iml When the vdW 
functional is included, the energy of the H-bond interac- 
tion increases by approximately 25 meV independently 
of the GGA used to describe the exchange interaction. 
Therefore, the weakening of the H-bond network by vdW 
interactions (reduction of the first and second coordina- 
tion peaks) is not associated to the weakening of the H- 
bond itself, but to the increase of new, favorable, non 
H-bonded configurations that compete with the H-bonds. 



B. Spatial Distribution Functions 

RDFs give angular-integrated information, but the an- 
gular distribution of molecules around a given one also 
contains very valuable information that can differentiate 
between similar RDFs. In Fig. [5] we plot the 0-0 spa- 
tial distributions functions (SDF)^- gooif, ^: 0) for three 
functionals. The polar angles 6, <j) are referred to a local 
coordinate set of the central molecule, with origin at the 
oxygen atom: x (direction of the HOH angle bisector), 
y (perpendicular to x, within the molecular plane) and 
z (normal to the molecular plane). We plot isosurfaces 
goo{T,9,(j)) = gc, restricted to spherical shells of thick- 
ness Sr = 0.2 A centered at three different distances r: i) 
at the first maximum of gooi^)^ with gc = 2, Figs. [5Ja- 
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FIG. 7: Total energy of the water dimer as a function of the 
intermolecular separation for two different molecular orienta- 
tions calculated for PBE (circles), revPBE (squares), DRSLL 
(diamonds) and DRSLL-PBE (triangles), (a): Non H-bonded 
configuration as shown in the inset (with partially facing O 
lone pairs from each molecule), (b): H-bonded configuration 
as shown in the inset. In both graphs the energies have been 
shifted to have the zero at the largest separation. 
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on the accepting molecule. 

The most important differences between the vdW and 
GGA functional occur beyond this tetrahedral shell, at 
the distance of the first minimum of .goo('') (Figs-IMd-f)). 
At this distance, the angular distribution for the PBE 
functional has small lobes in the tetrahedral directions 
and larger ones in the opposite directions, that we will 
refer as anti-tetrahedral;^ The small tetrahedral lobes 
correspond to the tail of the first coordination peak. The 
anti-tetrahedral lobes correspond to the directions of the 
interstitial sites, and it is fully dominant in the vdW-DFs 
(FigS- IHe-f)). Thus, although this anti-tetrahedral shell 
is already present with the GGA, due to entropic effects, 
its density largely increases in the vdW-DF due to the 
vdW attraction. Soper et al.-^ observed a similar effect 
in liquid water under pressure. Thus, the effect of non 
local vdW correlations in DFT water is similar to the 
effect of pressure in experimental water. 

The third shell shown in Figs. [5l^g-i), corresponds to 
the position of the second coordination shell, at r '^ 4.4 
A. This peak, formed by second neighbors, tetrahedrally 
coordinated to the first coordination shell, is depleted 
in vdW water in favor of the interstitial anti-tetrahedral 
shell. 

We can conclude from this analysis that the effect of 
vdW-DF correlation is to increase the density of the liq- 
uid by populating the interstitial sites that are at a dis- 
tance between the first and second coordination shells. 
These interstitial structures become more favorable be- 
cause of non-H-bonded vdW interactions, as shown in 
Fig.El^a), and they compete with, and destabilize, the H- 
bond network. They also increase the height of the first 
minimum in (?oo('')i that is too low with PBE. On the 
other hand, revPBE water is already relatively unstruc- 
tured, H-bonds being rather weak within this functional, 
as shown in Table IIIII Therefore the density increase 
induced by vdW interactions, combined with a weak H- 
bond network, results in a collapse of the structure, char- 
acteristic of water under pressure. 



C. Hydrogen Bond Network 



c); ii) at the first minimum of 5oo('')j with gc = 0.5, 
Figs. [SJd-f ) ; and in) at the second maximum of goo(''), 
with gc = 1, Figs. [5{g-i). To ease the visualization we 
show three different viewpoints (front, top and side) for 
each shell. 

Fig. E compares the 0-0 SDF of PBE, DRSLL, and 
DRSLL-PBE water at 1.0 g/cm^. We omit revPBE be- 
cause of its close similarity to PBE. The structure of the 
first maximum is almost identical for all the functionals 
considered, i.e, nearly tetrahedral. The density of ac- 
ceptor molecules (lobes in front of the H atoms) is much 
more localized than the density of donor molecules (lobes 
in the region of the oxygen lone pairs behind the oxygen 
atom) which is more disperse. This dispersion refiects 
the asymmetry of the H bonds, which are less directional 



In previous sections we did not explicitly analyze the 
hydrogen bond network. In water, each molecule tends to 
be surrounded by four others, donating two H-bonds and 
receiving two others in a nearly tetrahedral arrangement. 
This translates into a H-bond network with a majority of 
fourfold-coordinated molecules. However it also contains 
under-coordinated and over-coordinated molecules, and 
the population of under-coordinated molecules correlates 
with an increasing self diffusion constant. The average 
bond life-time of the H-bonds in liquid water is of the 
order of 1 psi^ 

One of the most striking differences between the GGA 
functionals and their vdW-DF counterparts is the large 
change observed in the first peak of the 0-0 RDF. In 
average, this implies both a reduced number of H-bonds 
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FIG. 8: Spatial distribution function(SDF) of liquid water, from 3 viewpoints (front,top and side) and for 3 different shell 
radii as indicated in the legend: (a) PBE, first coordination shell, (b) DRSLL, first coordination shell, (c) DRSLL-PBE, first 
coordination shell, (d) PBE, interstitial shell, (e) DRSLL, interstitial shell, (f) DRSLL-PBE, interstitial shell, (g) PBE, second 
coordination shell, (h) DRSLL, second coordination shell, and (i) DRSLL-PBE, second coordination shell. 



and longer H-bond lengths when vdW interactions are 
accounted for. This however does not correlate with the 
average H-bond energies and geometries presented in ta- 
ble mil The table shows that vdW interactions increase 
the H-bond energy by ~25 meV, both for revPBE and 
PBE. This is close to the vdW (non-H-bond) binding 
energy shown in Fig. El^a). Also, the right shift of the 
0-0 RDF peak (2.5%) is much larger than the increase 
in average 0-0 distances shown in table Hill (0.9%). To 
study the origin of these deviations we have analyzed the 
statistics and geometries of: (i) H bonds, defined by the 
molecules that contribute to the first, tetrahedral, shell 
shown in Fig. ISJa-c) and (ii) VdW bonds, or interstitial 
molecules in anti-tetrahedral coordination to the H-bond 
network as shown in Fig. [SKd-f). The total gooir) has 
been divided in three contributions: first-shell and sec- 
ond shells of H bonds (HBl and HB2) and a reminder 

Results are shown in Fig. M and Table |lV] for PBE 
and DRSLL-PBE, and in Supplementary Information 
for revPBE and DRSLL. VdW correlations produce 
slightly less H-bonds per molecule (3.78 vs 3.81) with 
a broader distribution of lengths. The second shell of 
H-bonds, represents ^ 60% of the second peak in the 
RDF. Again, there are slightly less second-H-bonded 
molecules with DRSLL-PBE (10.47 vs 10.59). These 
numbers are very close to the ideal 'mean-field' value 
Nhb2 = Nhbi{Nhbi - !)■ The bottom plot in Fig. M 
represents all the neighbors not included in the previ- 
ous two distributions. These include the third shell of 
H-bonds and the vdW-bonded molecules (these two are 
not mutually exclusive). We see an important 20% in- 
crease in DRSLL-PBE, up to a distance of 4.5 A, a region 



in which non-H-bonded interstitial molecules are domi- 
nant. This increase is about twice as large from revPBE 
to DRSLL (see Supplementary Information). Notice that 
the maximum of this distribution for DRSLL-PBE, at 3.5 
A, coincides approximately with the position of the vdW 
energy minimum of Fig. [3 and also with the first mini- 
mum of the total 5oo('')- We have also observed a pos- 
itive correlation between the number of non-H-bonded 
neighbors of a molecule and its average H-bond distances. 
This produces a 1.7% right shift of the average 0-0 dis- 
tance for the first H-bonded shell from PBE to DRSLL- 
PBE, due to the increase of vdW bonds. A more detailed 
analysis of this, and other related effects, will presented 
in a future contribution. 



TABLE IV: Mean values for H and vdW bonds at 1.0 g/cm^ 
and ~300 K. Nhbi and Nhb2'- average number of H-bonded 
molecules in the first and second coordination shells (integral 
of first two curves in Fig. |9]). NvdW- average number of non- 
H-bonded molecules (vdW-bonded and others) in a sphere of 
4.5 A radius, thbi'- average O-O distance for first H-bonded 
shell. 



Functional A'^;, 



Afj, 



N^dw Thbi (A) 



PBE 
DRSLL-PBE 



3.81 
3.78 



10.59 
10.47 



3.25 
3.89 



2.847 
2.895 



V. CONCLUSIONS 

We have performed DFT-based AIMD simulations of 
liquid water with different GGA and vdW functionals 
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FIG. 9: Decomposition of goo{r) into H and vdW bonds. 
Top, tetrahedral shell of first H-bonded neighbors, shown in 
Fig.|8|a-c). Center, second H-bonded shell shown in Fig.[8|g- 
h). Bottom, remaining molecules which include vdW-bonded, 
interstitial and third H-bonded molecules. The maximum 
of the lower DRSLL-PBE curve, at r=3.5A, is due to the 
molecules in anti-tetrahedral coordination to the H-bond net- 
work, as shown in Fig. [Sjd-f ) . Simulations in this figure are 
for liquid water at 1.0 g/cm^ and ~300 K. 



and different densities, comparing the structural and dif- 
fusive properties. The main conclusions are: 

• Liquid water under standard atmospheric pressure, 
as obtained using standard GGAs, is much less 
dense than experiments: PBE 12% less, BLYP 19% 
less and revPBE 31% less. Accordingly, the aver- 
age calculated pressure at the experimental density 
is as high as 1 GPa. 

• The net effect of vdW non-local correlations is to 
reduce the pressure and to increase the equilibrium 
density, which is just 2% larger than experiments 
with DRSLL and 13% higher with DRSLL-PBE. 

• With vdW interactions, non-H-bonded water-water 
configurations are introduced, due to the occu- 
pation of interstitial sites. These configurations 
have an anti-tetrahedral orientation, opposite to El- 
bonds, and they contribute to increase the height of 



the first minimum in (?oo(^)- While the interstitial 
sites are partially occupied due to entropy in GGA 
water, their population increases greatly when the 
vdW attraction is added. The use of angular- 
resolved spatial distribution functions is necessary 
to understand the overall structural changes. 

• These anti-tetrahedral structures are the key factor 
to increase the self diffusion of AIMD water. The 
diffusivity obtained with DRSLL is very close to 
experiment at room temperature, without the need 
of temperature rescaling. 

• DRSLL water is under-structured when compared 
to experiments, with a significant collapse of the 
second coordination shell induced by the increased 
density. This is related to the use of revPBE 
as local exchange, which produces too weak H- 
bonds. The combination of non local vdW correla- 
tion, as in DRSLL, combined with PBE exchange, 
produces a better liquid, with structural and dy- 
namical properties not far from experiments, with- 
out the need of any temperature rescaling. How- 
ever PBE H-bonds are too strong and DRSLL-PBE 
causes an over correction of the density with respect 
to PBE. 

• A better treatment of the exchange interaction 
should improve the H-bond description and, to- 
gether with vdW correlations, might end up provid- 
ing the final overall agreement with experimental 
results in terms of structure, density, and diffusivity 
of liquid water. A modification of the vdW density 
functional has been proposed very recently,— and 
a study of this and other flavors of vdW, not con- 
sidered here, are under study and will be published 
in a future work. 
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